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ABSTRACT 

The main accretion phase of star formation is investigated in clouds with 
different metallicities in the range of 0 ^ ^ Zq, resolving the protostel- 

lar radius. Starting from a near-equilibrium prestellar cloud, we calculate the 
cloud evolution up to ^ 100 yr after the first protostar formation. The star for¬ 
mation process considerably differs between clouds with lower {Z ^ 10“"^ Zq) 
and higher {Z > 10“^ Zq) metallicities. Fragmentation frequently occurs and 
many protostars appear without forming a stable circumstellar disc in lower- 
metallicity clouds. In these clouds, although protostars mutually interact and 
some are ejected from the cloud centre, many remain as a small stellar cluster. 

In contrast, higher-metallicity clouds produce a single protostar surrounded 
by a nearly stable rotation-supported disc. In these clouds, although frag¬ 
mentation occasionally occurs in the disc, the fragments migrate inwards and 
finally fall onto the central protostar. The difference in cloud evolution is due 
to different thermal evolutions and mass accretion rates. The thermal evo¬ 
lution of the cloud determines the emergence and lifetime of the first core. 
The first core develops prior to the protostar formation in higher-metallicity 
clouds, whereas no (obvious) first core appears in lower-metallicity clouds. 
The first core evolves into a circumstellar disc with a spiral pattern, which 
effectively transfers the angular momentum outwards and suppresses frequent 
fragmentation. In lower-metallicity clouds, the higher mass accretion rate in¬ 
creases the disc surface density within a very short time, rendering the disc 
unstable to self-gravity and inducing vigorous fragmentation. 

Key words: accretion, accretion disks—binaries: general—cosmology: theory— 
early universe—stars: formation—stars: protostars 
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1 INTRODUCTION 


Stars control the dynamical and chemical evolution of the universe. The typical stellar mass, 
stellar feedback and hnal fate of stars are closely related to the evolution of the universe, 
and the latter two factors are determined by the stellar mass at the moment of stellar 
birth. Thus, it is very important to understand the star formation processes throughout the 
history of the universe. Recent studies have shown that the star formation process in pri¬ 


mordial gas clouds considerab 


y differs from that i n present-day clouds and in clouds formed 


in the middle-aged universe ( 

Bromm et al. 

2001; 

Smith & Sigurdsson 

2007; 

JaDDsen et al. 

2007 

; 

Vlachida 

2008; 

Clark et al. 

2008 

JauDsen et al. 

200£ 

)l: Smith et al. 

2009; 

Machida et al. 

2009a 

DoDcke et al. 

2011, 

2013; 

Safranek-Shrader et al. 

2014 

). A major difference between 


primordial and present-day clouds is the abundance of metals and dust. Primordial gas 
clouds contain neither metals nor dust, while present-day clouds contain both. The abun- 


ance of metals and dust inf 


(lOmukai 


2000 


Omukai et al 


uences the t hermal evolution of collapsing star-forming clouds 


2005 


201011 . 


A collapsing gas cloud with a primordial composition cannot cool to below ~ 200 K, 
because, except for molecular hydrogen, there is no effective coolant at temperatures as low 
as T ^ 200 K. On the other hand, a present-day gas cloud can cool to ~ lOK because 
cooling by the metal and dust is effective even at T ^ 200 K. Numerical simulations have 


10^— 10° cm in the early universe (e.g. 


mi-halos) are formed at a num 

Der c 

ensitv of n-R ~ 

Abel et al. 

2002 

Bromm et al. 

2002 

Yoshida et al. 


20061) . while observations of nearby star-forming regions ha ve shown that molecula r cloud 


cores have a number density of lO'^ — 10^ cm ^ (e.g. iBergin fc Tafallal 120071) . Thus, 


star-forming clouds have almost the same density of n ~ 10^ — 10^ cm“^ in both the early 
and present-day universes. How ever, these clouds ha ve different temperatures with different 
abundances of metals and dust fjOmnkai et al.ll2005l ). With a cloud density of n lO'^cm 
primordial and present-day clouds have temperatures of ~ 200 K and ~ 10 K, respectively. 
If an equilibrium state is assumed to be realized in a star-forming cloud, the cloud mass 
depends on the gas temperature (or metallicity) under the assumption of spherical symmetry. 
Reflecting the difference in the gas temperature in clouds, star-forming clouds at each age 
should have different masses. For example, primordial star-forming clouds have a mass of 
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Different Metallicities 3 


^ 10^ flBromm et al.ll2002l) . while present-day clouds (molecular cloud cores) are ~ 1 Mq 
f Andre et al. 2010 1. 

After protostar formation in a collapsing cloud, the mass accretion rate onto the proto¬ 


star also depends on the gas temperature of the star-formi ng cloud. 


accretion rate is expected to be proportional to ~ c^/G (jShul Il977 : 


leoretically, t he mass 


Larson 


20031), where 


Cs is the sound speed of the prestellar cloud. Thus, the mass accretion rate should be ap¬ 
proximately 100 (200K/10K)^'®) times larger in a primordial cloud than in a present-day 

cloud. T herefore, primordial stars are considered to be more massive than present-day stars 


(see also 


Omukai Sz Palla 


2001 


20031) . However, this rough estimate ignores fragmentation 


in the collapsing cloud. A single massive star forms without fragmentation, whereas vigorous 
fragmentation leads to a multiple stellar system comprising less massive stars. 

We cannot directly observe the formation sites of stars even in present-day clouds be¬ 
cause they are embedded in dense gas envelopes. Thus, the star formation process must 
be clarihed by a theoretical approach and numerical simulations. Several numerical simu¬ 
lations have investigated the evolution of primordial and present-day star-forming clouds. 
However, simulating a collapsing cloud over sufficiently long durations to resolve the proto¬ 
star is very difficult, because the numerical timestep shortens as the gas density increases. 
Since the protostar has a very high density, cloud evolution around a protostar requires an 
extremely short timestep precluding long-term integration. The sink method is used to over¬ 
come this problem, in which the high- density gas is assumed to become a protostar and is 


repla ced by sink particles or sink cells flBate et al.lll995 : 


Krnmholz et al. 


2004 


Hubber et al. 


20131 ). At the expense of spatial resolution inside and near the sink accretion radius, some 


stud ies have realized the long-term Simula 


ties fBromm et al. 


ion of star-forming clouds with different metallici- 


20Q1; 

Clark et al. 

2008; 

Dopcke et al. 

2011, 

20ls3; 

Safranek-Shrader et al. 


20141) . Although these studies have provided insight into large-scale fragmentation, star for¬ 


mation efficiency and stellar mass distribution, the evolution of the circumstellar region is 
masked by the sink. In the protostar proximity, we can expect small-scale circumstellar disc 
formation and fragmentation, which may affect the large-scale cloud evolution. 

Some studies have inves tigated 


withou 

; emploving a 

sink 

2 nnb: 

Tomida et al. 

2013 


Bate 


1998 


Bate et al 


Machida et al. 


in present-dav star- 

brming clouds 

2997: 

Bate 

2911: 

Machida et al. 


20141 ). These studies could model cloud evolution 


for only ~l-100yr after the protostar formation, and demonstrated that fragmentation 
is rare in high-density gas regions. In addition, these studies indicate that a rotation- 
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supported disc ioims befor e the protostar formation, where the hrst core flLarson 


1969; 


Masunaea fc Inutsukal 120001) formed before the protostar formation evolves into a circum- 


stellar disc. Thus, at the protostar formation epoch, the protostar is already surrounded 
by a large-scale rotation-supported disc with spiral structures, which effectively transfer 
angular momentum and stabilize the disc against gravitational instability and subsequent 
fragmentation. 


On the other hand, the evolution of primordial star-formin g clouds on protostar-resolving 


scales 


291811 . 


las been scarcely i nvestigated in models without a sink (iGreif et ah 


2012 


Machida fc Doi 


Greif et ah 


(1201211 showed that frag mentation is common around initially formed pri¬ 
mordial protostars. iMachida fc Poll (1201311 reported that no stable circumstellar disc is es¬ 
tablished around a primordial protostar because vigorous fragmentation and the orbital mo¬ 
tion of the fragments signihcantly disturb the circumstellar environment. They also showed 
that when a sink is introduced into a primordial cloud, the fragmentation scale is (arti- 
hcially) determined by the sink accretion radius. Furthermore, in primordial star-forming 


hrst core does not appear in the primordial collapsing gas ( 

Omukai 

O 

O 

O 

CM 

Omukai et ah 

2005 

Yoshida et ah 

2QQb; 

Greif et ah 

2Q12, 

2 QlsT 

Machida & Doi 

2013 

)• 


Previous studies without a sink reveal signihcant differences between the primordial and 
present-day star formation processes. Primordial clouds frequently fragment without form¬ 
ing a stable circumstellar disc, whereas the stable rotation-supported disc prevents frequent 
fragmentation in present-day clouds. Thus, it is considered that the star formation process 
depends upon the cloud metallicity. The large-scale cloud evolution after the protostar for¬ 
mation at different metallicities has been widely investigated with a sink, but has not been 
investigated without a sink. Although a sink may be required in long-term cloud evolution 
calculations, we should hrst study the cloud evolution by resolving the protostar to investi¬ 


gate the effect of the small-scale structure around the protostar on the large-sca. 


In this study, we adopt a protostellar model similar to that of 


Machida fc Doi 


e struc ture. 
J 2 OI 3 I 1 and 


investigate the star formation process from the prestellar cloud stage until approximately 
100 yr after the protostar formation, varying the cloud metallicities within the range of 
0 ^ Z ^ Zq. 

The rest of the paper is structured as follows. Our model framework and numerical 
methods are described in §2 and the protostellar model is explained in §3. The numerical 
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Different Metallicities 5 

results are presented in §4, while §5 is devoted to the fragmentation condition and effects of 
the initial condition. The results are summarized in 56. 


2 MODEL AND NUMERICAL METHOD 


We calculate the evolution of eight clouds with different metallicities in the range of 0 ^ 


Z ^ Zq. The initial clouds have the critical Bonner-Ebert (BE) d ensity proh 


is ch aracterized by the central density and isothermal temperature flEbert 


1955 


e , which 


Bon nor 


19561) . To construct the BE profile, we set the central density as n = lO'^cm”^ in all mod¬ 
els and adopt different isothermal temperatures that are estimated by a one-zone calcula- 
tion (see below). To promote contraction, we enhance the cloud density by a factor of 1.8 
flMatsumoto fc Hanawall2003l) . For each initial cloud, the mass and radius are different, while 
the initial cloud stability (the ratio of thermal to gravitational energy), which affects the 
fragmentation condition, is the same (for details, see ^5.2p . The cloud metallicity Z, initial 
(central) temperature T^, cloud mass M^i and radius Rd in each model are listed in Tabled! 

For all models, a rigid rotation of Hq = 3 x 10“^^ s“^ is adopted. Each initial cloud has 
the same ratio of rotational to gravitational energy of /So = 10“^. Although the rotational 
energy adopted in this study lie s within the observa tional constraints of nearby star-forming 

I 2 OO 2 ) , it may be somew hat smaller than that 


Caselli et ah 


regions (e.g. 10““^ < /Sq < 0.07; 
estimated by recent cosmological simulations (e.g. 
value of /So for all models to limit calcu 
influence the fragmentation process (e.g. 


Hirano et ah 


20141) . We adopt the same 


ation models, while the initial rotation rate may 


Machida et ah 


20081 ). In addition, the initial cloud 


mass may also be important for investigating the star formation process. The effects of the 
initial rotational energy and cloud mass are discussed in 55.21 


T he cloud eyolutiqn is com puted by the nested grid method (for details, see 


200 


Machida et ah 


2005a 


Machida et ah 


b|), which constructs nested rectangular grids with the same 


number of cells (i, j, k) = (256, 256, 32) or (128, 128, 16). We impose a mirror symmetry 
with respect to the z = 0 plane. Initially, we prepare the hfth grid level and immerse the 
initial cloud in the hrst grid level, whose box size is twice the initial cloud radius 2Rc\. A 
hxed boundary condition is imposed on the surface in the hrst grid level. As the calculation 
proceeds , hner grids are dyna mically generated, which resolve the Jeans length in at least 


16 cells (ITruelove et ah 1119971 ). Although the grid sizes and cell widths differ between the 


models even at the same grid level /, we can appropriately calculate the cloud evolution 
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and fragmentation process because the Truelove condition is always fulfilled (except for the 
escaped fragments, see § 1 ]). 

Our numerical method solves the following hydrodynamics equation with self-gravity: 
dp 


dt 


+ V ■ {pv) = 0, 


dv 


p— + p{v ■ V)v = -VP - pV0, 

= AnGp, 


( 1 ) 

( 2 ) 

(3) 


where p, v, P and (j) denote the density, velocity, pressure and gravitational potential, re- 
spectivel y. For gas pres s ure, w e adopt a barotropic relation from the one-zone calculation 


given by 


Omukaietah 


(12005^ . Their one-zone model calculates the thermal evolution of 


collapsing gas clouds with different metallicities (0 ^ ^ ^o)- In Figure [H the thermal 

evolution of a collapsing gas cloud at different metallicities is plotted against th e number 


density. An almost identica l barotropic relation was used in our previous studies flMachida 


2008 


Machida et ah 


2009al) . In the present study, we adopt a slightly different thermal evo¬ 


lution. Specihcally, we impose a lower temperature limit of 10 K to model a more realistic 
scenario. In previous studies, t he gas in clouds with Z = cools to ^ 3 K in the range 
of lO^cm”^ ^ n ^ 10 ^°cm“^ f Machida 2008 : 


Machida et ah 


2009al) . Note also that since 


our present study focuses on high-density gas and the early star formation process, this 
difference does not signihcantly affect the results. 

For each metallicity, two evolutionary tracks are plotted in Figured! One track is given 
by the one-zone calculation, which yields a smoothly connected evolution from low- to high- 
density gas. In the other track, the protostellar model is added to the original one-zone result. 
The protostellar model visibly affects the results only at high density (n 3x10^^ cm“^). The 
polytropic index 7 of the original one-zone calculation is ~ 1.1 in the range of lO^^cm”^ ^ 
n ^ 10^*^ cm“^, whereas that of the protostellar model is ~ 4 in the high-density gas region. 
The protostellar model is described in the next section. 


3 PROTOSTELLAR MODEL 

In many previous studies, the sink method is introduced to investigate long-term cloud 
evolution, in which the high-density gas region is masked by sink cells or sink particles and 
is not spatially resolved. However, as described in [HI a sink tends to mislead the result 
because, for example, the fragmentation scale is artihcially determined by the sink accretion 
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radius flMachida fc Doill2013() . To avoid artificial effects caused by a sink and to investigate 


the early star formation process in a 
model, the same as that adopted in 


ligh-density gas r egion , we here use a simple protostellar 


Machida fc Doil (120131) . where the thermal evolution or 


the polytropic index is adjusted in the high-density gas region and the protostellar radius is 
related to the protostellar mass. 


To connect the proto stellar mass an c 


relation demonstrated by 


Omukaietah 


rad ius, we refer to the protostellar mass-radius 


(120101) ■ who used their one-zone calculation to es¬ 


timate the mass accretion rate onto the protostar in clouds with different metallicities. In 
our study, we first check the long-term mass accretion rate by setting up non-rot at ional 
equivalents of the initial clouds listed in Table [T] and by calculating them with a sink. We 
set the sink accretion radius to race = lAU and the threshold density to nthr = 10^®cm“^. 
Then, we calculate the cloud evolution for ^ 10^ — 10^ yr after the protostar formation, in 
which we define the protostar formation epoch at which t he maximum cloud de nsity reaches 
the threshold density (for details of the sink method, see lMachida fc Doill2013l) . 


The mass accretion rate for different metallicities is plotted against the protostellar 
mass in Figure O In clouds with higher metallicity {Z > 10“^ Zq), the mass accretion rate is 
small, and hence, we can calculate the cloud evolution up to a protostellar mass of ~ 0 . IMq. 


Although we only calculate the very early phase of star 
in Fi gure [2] approximate those reported in 


Omukai et al 


brmat ion, the mass accretion rates 


( 2ninli 


(see also 


Tanaka Sz Omukai 


2014J) . Since the accretion history mainly governs the protostellar evolution, we expect that. 


in our calcula t ion, t he protostar evolves similarly to the evolutionary process described in 


Omukai et al. 


(120101 ). Thus , we can relate the p rotostellar radius to the protostellar mass 


by referring to Figure 17 in 


Omukai et al. 


(I2010h . 


Next, to construct the protostell ar model, we ca. 


cloud without a sink, as described in iMachida fc Doil (120131 ). parameterizing the polytropic 


c ulate the evolution of a non-rotating 


index 7pg and the protostellar density Ups at which the thermal evolution differs from the 
one-zone result. For each model listed in Table [H we calculate the protostellar evolution 
of 5-20 clouds, varying Ups and 7ps, and determine the most pl ausible parame t ers of the 


protostellar model, with which the mass-radius relation given by 


Omukai et al. 


( 201ol) 


are 


best reproduced. Figure [T] plots the thermal evolution of the protostellar model (red lines) 
at each metallicity. As shown in the figure, the protostellar densities are in the range of 
3 X lO^^cm”^ ^ Ups ^ 3 X lO^Wm”^. The polytropic index of 7 ps = 4 is adopted for all 
models. The protostellar density in each model is listed in Table [TJ With a harder equation of 
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state, long-term calculation becomes possible because the density increase almost stops inside 
the protostar (or inside the shock surface), alleviating the Truelove condition. Although we 
cannot resolve the inner structure of the protostar, we can resolve the circumstellar region 
just outside the protostar. Note that we do not resolve the high-density gas region (n Ups) 
with a stiff EOS, especially for lower-metallicity models {Z ^ 10 Although the mass 

accretion rate is mainly deter mined by the large sca le structure without a persistent Erst 
core formation in such models (iMachida et al.ll2009al ). a stiff EOS may somewhat influence 
the mass accretion rate and fragmentation process. 

Figure |3] plots the protostellar radius against the protostellar mass for different metallic- 
ities, in which the cloud evolution is calculated with our protostellar model (Fig. [1]) and the 
protostar is dehned as the object surrounded by the shock surface in a high-density gas re¬ 
gion of n > 3 X 10^^ c m~^. In the hgure, th e protostellar radii for Z = 10“®, 10““^, 10“^, 10“^ 


and 1 ^0, taken from 


Omukai et ah 


(120101) ■ are also plotted by thin li nes. The hgure show s 


that our protostellar model reproduces well the mass-radius relation in 


Omukai et al. 


( 120101 ) . 


The hgure indicates that protostars formed in clouds with Z = 0-10“® Zq have a radius of 
30 Rq ^ Rps ^ 200 Rq when the protostar has a mass of Mpg < 10 Mq. On the other hand, 
the protostellar radii are in the range of < 30 Rq for models with Z ^ 10 The proto- 


Omukai et al 


stellar mass-radius relations in Figure |3] are in good agreement with those in 
(120101) . The protostellar radius increases as Rps oc during the adiabatic accretion 


phase. This phase lasts unt i 


(Hosokawa fc Omukai! 


20091: 


the protostar ente r s the Kelvin-Helmholtz contraction phase 


Omukai et al 


2010 ; 

Smith et al. 

2011 , 

2012 ) 


20121 ) . As indicated in 


Figure [HI the protostellar radius in each model is roughly proportional to Rps oc Mp/^. Thus, 
the protostellar model adopted in this study roughly reproduces the protostellar evolution 
(or the mass-radius relation) reported in previous studies. With this protostellar model (or 
EOS), we calculate the evolution of clouds with diherent metallicities. To check the pro¬ 
tostellar evolution, the detailed evolution of non-rotating clouds is described in Appendix 


4 RESULTS 

To investigate the ehects of metallicity on the star formation process, we prepared prestel- 
lar clouds, varying the metallicity in the range of 0 ^ Z ^ ^0, and calculated them for 
~ 100 yr after the first protostar formation (or ~ 10® yr after the initial cloud begins to 

© 0000 RAS, MNRAS 000, 000-000 

































Different Metallicities 9 

collapse). During the calculations, in any model, fragmentation occurred in the region of 
r < 20 AU and we could cover the fragmentation region by the hnest grid. However, some 
fragments escaped from the central region by mutual gravitational interaction between frag¬ 
ments, and fragments that had moved far off-centre could not be resolved to sufficiently 
high spatial resolutions (see §2]). Thus, we dehned the fragments separated by more than 
20 AU as ‘escaped fragments’ and ignored their evolution in subsequent analysis. Note that 
although the escaped fragments were calculated with relatively low spatial resolution, they 
did not signihcantly affect the cloud dynamics because their mass 0.1 Mq) was negligible 
relative to the total cloud mass (7 — 2000 Mq). 

In the higher-metallicity models, we have to resolve a smaller structure or a higher density 
region because a protostar formed in a high-metallicity cloud has a smaller size and a higher 
protostellar density Ups, as described in l|3l To spatially resolve such small protostars, we 
require a hner grid with smaller cell width around the cloud centre. Thus, long-term cloud 
evolution is limited in these models because a hner resolution requires a shorter time step 
and much longer CPU time. However, for the high-metallicity models of Z ^ Zq, we 
conhrmed that fragmentation is rare and that high-density gas and the protostar(s) occupy 
the limited region around the cloud centre, as described in the following subsections. In such 
models, a large number of cells are not required in each nested grid to maintain the Truelove 
condition. In the Z = 10“^, 10“^ and 1^0 models, we calculated the cloud evolution with 
cell numbers (i, j, k) = (256, 256, 16) for ~ 30 yr after the hrst protostar formation, and 
conhrmed that fragmentation occurred only in the region near the hrst formed protostar (or 
in the circumstellar disc around the protostar). To reduce the CPU time and enable long¬ 
term calculation, we then calculated the cloud evolution with cell numbers {i, j, k) = (128, 
128, 8), which is sufficient to maintain the Truelove condition. In this setup, we calculated 
the cloud evolution for over 100 yr after the hrst protostar formation in the Z = 10“^ 
and 10“^ Zq models, in addition to the Z = 0, 10“®, 10“®, 10“^ and 10~^ Zq models. The 
computational time consumed by each model was between four and six months of wall clock 
time. In the Z = \ model, we only calculated the cloud evolution for 55 yr after the hrst 
protostar formation because a considerably high spatial resolution and considerably longer 
calculation time were required. 

In this section, we hrst describe the evolution of each cloud in the main accretion phase 
during which the protostar signihcantly increases its mass by mass accretion from the in¬ 
falling envelope. Then, we compare the circumstellar structure, number of fragments and 
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protostellar mass between the different models. Although we calculated the collapsing cloud 
from a prestellar stage with a central density of 10^ cm“^, our main focus is the cloud evolu¬ 
tion after the hrst protostar formation (i.e. the main accretion phase). The evolution of the 
gas-collapsing ph ase prior to the protostar formation has already been investigated in our 


previous studies (IMachida et ah 


2006 


2008 


Machida 


2008; 


Machida et ah 


2008c . 


2009a 


In the following subsections, we calculated the cloud evolution em with our protostellar 
model without a sink. 


4.1 Z = 0 and 10 ® Zq Models 


In this subsection, we briefly describe the cloud evolution in the Z = 0 and 10’ models. 
The cloud evolution in the Z = 10 mo del is almost identical to that in the primordial 
cloud, which is detailed in IMachida fc Poll (l2013f) . Figure 0] shows snapshots of the Z = 0 
(left) and 10~^ Zq (right) models at tps ~ 100 yr, where tps is the elapsed time after the 
hrst protostar formation and tps = 0 is dehned as the epoch at which the maximum density 
reaches n^ax = ''^ps (see Tabled]). In the hgure, the density (upper panels) and temperature 
(lower panels) distributions are plotted, in which the temperature distribution is constructed 
using the density distri buti on and barotrop i c equ ation of state. Resolving the protostellar 


scale. 


Greif et al. 


(120121) and IMachida Sz Doil fl2013f) showed that collapsing primordial clouds 


frequently fragment and that the resulting fragments (or protostars) tend to cluster around 
the centre of the collapsing primordial cloud. We refer to a fragment as a protostar when 
its maximum density in the fragment exceeds the protostellar density, i.e. when Umax > n 


pS5 


and as a clump when Umax < ’^ps- The upper panels 


simi 


2011 


a r to t h ose of previous studies (jClark et al. 


2012 : 


Machida fc Doi 


2008 


of Figu r e IH show results t 


2011a b 


Smith et al. 


lat are very 


2011 


Greif et al. 


20131) . Six fragments are clearly observed in the left panel of 


Figure m the most massive protostar has a mass of ~ 2Mq with a radius of ~ IGOR©. In 
the Z = 10“® Zq model, fragmentation frequently occurs after the hrst protostar formation 
and seven fragments appear, as seen in the Figured] right panel, in which the most massive 
protostar has a mass of 1.8 Mq with a radius of ~ 140 Rq. 

In both upper panels, the most massive protostar is surrounded by a hlamentary struc¬ 
ture. In addition, no stable rotation-supported disc appears after the hrst protostar for¬ 
mation in either model. Although a circumstellar disc transiently appears around massive 
protostars, it immediately fragments into less massive protostars. Fragmentation continues 
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in the central region of the clond and over 10 to 20 fragments are observed at most by the 
end of the calcnlation. The clond evolntion dnring the main accretion phase is essentially 
the same in the Z = 0 and 10 models because both models have almost identical ther¬ 
mal histories (Fig. [T]), mass accretion rates (Fig. [2]) and protostellar mass-radius relations 

(Fig. ED. 

The temperature distribution basically traces the density distribution. However, we can 
see an interesting feature around some fragments in the left panel {Z = 10~^ Zq model): 
a low-temperature region appears around or inside some fragments. The low-temperature 
region corresponds to a temporary temperature drop at n ~ 10^® cm“^ (see Z = 0 evolution 
track of Figure ED. This indicates that fragmentation tends to occur in a relatively low- 
temperature region. 


4.2 Z = 10-5 Zq Model 

Figure El shows the time sequence of the density distribution on the equatorial plane around 
the centre of the collapsing cloud in the Z = IQ-^ Zq model. The cloud evolution in this 
model is qualitatively similar to those in the Z = 0 and 10“® Zq models; however a slight 
quantitative difference is evident. As shown in Figure El fragmentation occurs just before 
the maximum density reaches rips (Fig. Eh), and two fragments appear (Fig. Eh)- Further 
fragmentation then occurs around the hrst formed fragments that have the maximum mass 
and size (Fig. Eh and d) among fragments. For tps 50 yr, a hlamentary structure on a 
scale of 10 — 30 AU develops and fragments over time. As seen in Figures Eh-5', more than 
10 fragments appear within the region of r < 30 AU, where r is the radius from the cloud 
centre. 

Less massive protostars tend to fall onto more massive protostars. Thus, in addition to 
gas accretion, the massive protostar acquires its mass by mergers with fragments. Moreover, 
less massive protostars are ejected from the cloud centre by mutual gravitational interaction 
between protostars. At the end of the calculation, 11 protostars had been ejected from 
the cloud centre in the Z = IQ-^ Zq model. Throughout the remainder of the evolution, 
these protostars moved nearly radially outwards in the region of r > 20 AU. The number 
of fragments increased up to tps ~ 90 yr (Fig. E]f-^) ^md temporarily decreased thereafter 
(Fig. ED)- Fragmentation is ineffective around ma ssive protostar s that can stabilize their 
circumstellar gas against gravitational instability flToomrel 1 19641) . On the other hand, the 
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fragmentation region extends over time because gas with larger angular momentum falls 
later. Since we calculated the cloud evolution for only ~ 100 yr after the hrst protostar 
formation, we cannot determine whether vigorous fragmentation continues or whether the 
system stabilizes. However, we can conclude that fragmentation is a frequent event and that 
many protostars appear in the very early phase of star formation in the Z ^ 10“® Zq models. 

Figure E] left panel shows the temperature distribution on the equatorial plane at = 
60.6 yr and indicates that the central region of the collapsing cloud (r 20 — 40 AU) has 
a lower temperature than the outer envelope (r 20 — 40 AU). The inner low-temperature 


region corresponds to a temperature decrease in the range c 

)f 10^^ cm ^ ^ 10^® cm ^ of 

Figure [H which is caused by dust cooling ( 

Omukai et a! 

. 200 

5). A relatively low temperature 

can cause frequent fragmentation ( 

Tanaka & Omukai ^ 

>014) 



4.3 Z = 10 ^ Zq Model 

In the Z = 10“^ Zq model, fragmentation continues around the cloud centre until the end of 
the calculation, but the number of fragments is smaller than that in the Z ^ 10“® Zq models. 
Figure [7] shows the time sequence of the Z = 10“^ Zq model. The hrst fragments appear 
just before the protostar formation (Fig. [7]a). These fragments merge into a single protostar 
enclosed by a disc or spiral structure (Fig. ^). An additional two fragments appear in the 
disc, as seen in Figure [7]c. A hlamentary structure then develops, with further fragmentation 
within or near the hlament (Figs. Wd-g)- Some of smaller protostars fall onto more massive 
protostars, and some are ejected from the cloud centre. Only two protostars remain in the 
region of r < 5AU at tpg ~ 70 yr (Fig. ^). Subsequently, in the evolutionary period of 
fps ~ 85 — 100 yr, only one to three protostars and a few clumps (maximum density of 
’^max < ’^ps) remain. However, for ^ 100 yr, fragmentation is again activated and over 
hve protostars appear in the region of r < 20 AU (Fig. O). Figure [6] right panel shows the 
temperature distribution at t-pg = 63.2 yr and indicates that the fragmentation region is 
embedded in the lower-temperature region caused by dust cooling, as seen in Figured left 
panel. Comparing Figures 0] and E] with Figure [TJ the number of fragments is substantially 
smaller in the Z = 10“^ Zq model than in the Z < 10' models. 
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4.4 Z = 10 ^ Zq Model 


Unlike the Z ^ 10“'^ Zq models, fragmentation rarely occnrs in Z > 10“^ Zq models. In 
the Z = 10“^ Zq model, fragmentation occurs in the outer disc region and any fragment 
falls onto the central protostar before it can evolve into separate protostars (that is, before 
its maximum density reaches the protostellar density Upg). In this model, a single protostar 
(and a few clumps) remains at the end of the calculation. 

Figure ID shows the time sequence of the Z = 10’ model. As observed in the Z = 
10“^ Zq model (Fig.[7j), the hrst fragmentation in this model occurs just before the maximum 
cloud density reaches the protostellar density Upg, and two fragments appear (Fig. [H]a). 
Both fragments evolve into protostars that orbit each other (Figs. [8]& and c) before merging 
into a single protostar at fpg ~ 16 yr. A circumstellar structure and spiral pattern then 
develop around the protostar. After the merger, no further fragmentation occurs until fpg ~ 
70 yr (Figs. [H]d-/). Thereafter, some fragmentation occurs in the outer disc region (Fig. [Hb'). 
While such fragments (or protostars) are orbiting around the central massive protostar, they 
gradually migrate inwards under gravitational interaction with the circumstellar disc and 
hnally fall onto the central protostar. One of these fragments evolves into a protostar before 


falling onto the central protostar, whereas the others fall before 

protostellar density. Inward migration and falling onto the central 

their density reaches the 

protostar are considered to 

tvpifv the present-dav star formation processes ( 

Vorobvov & Basu 

2006 

2010 

Machida et al. 


2011 a 


Tsnkamoto fc Machida 


2011 


Tsnkamoto et al 


20131 ). The cloud evolution after the 


hrst protostar formation considerably differs between the Z = 10 ^ Zq and Z ^ 10 Z^ 


o 


models; the former is more representative of present-day star formation processes. 


4.5 Z ^ 10 ^ Zq Models 

The cloud evolutions and circumstellar environments in the Z = 10“^, 10“^ and 1^0 models 
are similar to those in the Z = 10’ model. Figure [9] shows the time sequence of the 
Z = 10“^ Zq model. Also for this model, we hnd that fragments appear just before the 
protostar formation (Fig. [9]a) and later merge into a single protostar (Fig.[9]5). The merger is 
followed by the formation of a stable circumstellar disc with no clear spiral structure (Figs. [9]c 
and d). The inner edge of this disc fragments at fps — 50 yr (Fig. [9]e), but the fragment 
immediately falls onto the protostar and disappears (Fig. [9](). Although fragmentation in the 
circumstellar disc continues up to tpg ~ 70 yr (Fig. [9]g), the fragments always move inwards 
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and merge with the central protostar. At the end of the calculation, a single protostar 
remains, enclosed by a stable rotation-supported disc (Fig. IHli). 

Figure dn] shows the time sequence of the Z = 10“^ Zq model. As observed in other mod¬ 
els, the hrst fragmentation occurs just before the hrst protostar formation (Fig. [Tllhh The 
fragments orbit each other for a while, then they merge into a single protostar. Although 
fragmentation occurs again near the protostar just after the merger (Fig. [TObl. the fragment 
immediately merges into the protostar fFig.fTOldh Note that in reality, fragmentation may not 


occur in the very proximity of the protostar because the protostef 


nored in this study, can hea t the disk and suppress fragmentation flStamatellos fc Whitworth 


ar luminosity, which is ig- 


2009aUb : 


Lomax et ah 


20141) . Thereafter, a spiral structure develops without further fragmen¬ 


tation, as shown in Figures flOle-i At the end of the calculation, a single protostar is enclosed 
by a rotation-supported disc of radius of ~ 5AU. 

Figure [TT] shows the density distribution in the Z = IZo model at tps = 45.1yr. As 
previously mentioned, we calculated the cloud evolution for only 55 yr after the protostar 
formation in this model. The hgure shows a single protostar enclosed by a circumstellar disc. 
The spiral pattern develops in the circumstellar disc, while no fragmentation occurs by the 
end of the calculation. 


4.6 Number of Fragments 

As mentioned in previous subsections, the frequencies of fragmentation and the numbers 
of fragments considerably differ between the models. Lower-metallicity clouds are prone 
to frequent and numerous fragmentation events, whereas higher-metallicity clouds tend to 
establish a single protostar enclosed by a stable rotation-supported disc. We quantitatively 
estimated the fragment numbers and protostellar masses among the models by the following 
procedure: (1) the maximum density Umax is searched for within the region of r < 20 AU; (2) 
when the maximum density exceeds the protostellar density Umax > ''^ps, the centre of the 
protostar (or protostar position) is assumed to be on the cell having the maximum density; 

(3) the protostellar mass is estimated by integrating the mass over n > 0.25 Ups within the 
region of Tp < 2 Aj, where rp and Aj denote the radial distance from the protostellar position 
and the Jeans length (derived from rips and the corresponding temperature), respectively; 

(4) the integrating region (n > 0.25 Ups within Rps < 2Aj) is masked and ignored in the 
next step of searching for the maximum density; (5) steps (l)-(4) are iterated until the 
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maximum density is less than n^s- The number of fragments and mass of each protostar are 
then determined. 

This procedure was implemented at different time steps for each model. Note that when 
counting the number of protostars, we omitted the protostars that escaped from the region 
of r < 20 AU, as mentioned above. Note also that the density and spatial ranges over which 
the protostellar mass was integrated, namely n > 0.25 Ups and rp < 2Aj, were decided by 
trial and error. In these trials, the number of fragments output by the numerical procedure 
was compared with that counted manually. Although this procedure roughly estimated the 
number of protostars and their masses, it has limited applicability. For example, very close 
binaries separated by < 2Aj are recognized as a single protostar. However, we consider that 
this procedure allows a qualitative, rather than a quantitative, discussion of the relation 
between the fragmentation process and cloud metallicity. 

Figure [12] plots the number of protostars in all models against the time elapsed after 
the hrst protostar formation. The hgure clearly indicates that the number of protostars 
decreases as the cloud metallicity increases. Although we calculated the cloud evolution 
for only ~ 100 yr after the hrst protostar formation, the effect of cloud metallicity on the 
protostar formation process is clear. 

No signihcant difference is apparent between the Z = 0 and 10“® Zq models; both models 
have at most 17 or 18 protostars around the centre of the collapsing cloud. More protostars 
appear in the Z = 10“^ Zq model than in the Z = 0 and 10“® Zq models. The Z = 
10“® Zq model yields at most 25 protostars. This difference may be attributed to the higher 
protostellar masses in the Z = 0 and 10 models than in the Z = 10 ^ Zq model f 1 14.71) . 
Massive stars can stabilize their circumstellar environment against gravitational instability 
(for details, see 1 15.1.2p . In addition, the temperature decrease by dust cooling can affect the 
number of fragments, as shown in Figures 01 and O Considerably fewer protostars develop in 
the Z = 10“^ Zq model than in the Z < 10-4 Zq models. Although 10 protostars appear in 
the Z = 10-4 Zq model at tps ~ 80 yr, most of them subsequently merge into a high-mass 
protostar, yielding a single star at tps ~ 90 yr. However, disc fragmentation yields several 
protostars around the most massive protostar for fps 100 yr. 

Fragmentation is rare in the Z > 10 models. Although two protostars are formed 
in the very early accretion phase (tpg ^ 10 yr) in the Z = 10-^, 10-^ and 10-^ Zq models, 
they merge into a single protostar. A rotation-supported circumstellar disc then develops 
around the merged protostar. Although the disc occasionally fragments, the fragments (or 
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protostars) eventually fall onto the central protostar. Consequently, these models always 
yield a single protostar at the end of the calculation. The Z = 1 Zq model also produces 
a single protostar with no fragmentation. In the Z > 10-"Zo models, the spiral pattern 
induced in the disc appears to facilitate outward transfer of angular momentum. In addition, 
fragmentation is likely suppressed by the lower mass accretion rate during the main accretion 
phase. 


4.7 Protostellar Mass 


Figure [13] plots the mass of the most massive protostar in each model against the time elapsed 
after the hrst protostar formation. Note that this hgure does not necessarily trace the same 
protostar because the most massive protostar among multiple protostars can vary from one 
epoch to another. Note also that since we integrated the high-density gas in a region fl l4.6p 
without a sink, the protostellar mass may be somewhat overestimated or underestimated. In 
this hgure, a sudden increase indicates merging of a less massive protostar, while a sudden 


decre ase can result from hssion or tidal gas stripping by another protostar flMachida fc Doi 


20131 ). The sudden increase and decrease are also caused by changes in the corresponding 


protostar. At tps ~ 100 yr, the most massive protostar is ~ 1 Mq and ~ 0.05 Mq in the 
Z ^ 10“^ Zq and Z > 10“^ Zq models, respectively. Thus, the mass is accreted at an 
estimated rate of ~ O.OlMQyr”^ in the Z ^ 10“^ Zq models and ~ 5 x 10“^MQyr“^ 
in the Z > 10-4 Zq models. The mass accretion rates in different classes (Z ^ 10’ 
and Z > 10-4 Zq) differ by a factor of ^ 10 and roughly agree with those in Figure [2l 
To understand the protostellar evolution, we need to calculate a further evolutionary stage 
because the protostellar mass is much less than the host cloud mass even at the end of the 
calculation. 

Figure [m shows the protostellar masses in the Z = 0, 10-®, 10-^, 10-4, 10-^ and 10-^ Zq 
models, in which the masses of the hrst, second, fourth, sixth, eighth and tenth most massive 
protostars are plotted against the time elapsed after the hrst protostar formation. The total 
mass of the protostars is also plotted by a broken line in each panel. Since the number of 
protostars is less than Apg = 10 in the Z > 10-4 Zq models, fewer than seven lines are plotted 
in the panels of Z = 10-^ and 10’ -^^0. This hgure indicates that the masses of the most and 
second-most massive protostars are very similar in the Z = 0, 10-® and 10’ models. In 
the Z = 10-4 Zq model, however, the most massive protostar is substantially more massive 
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than the others. Moreover, in the Z = 10“^ and 10“^ Zq models, the secondary protostar 
is considerably less massive than the most massive protostar for tps 60 yr. Note that in 
these models, the masses of the most massive and secondary protostars are comparable only 
when fragmentation occurs near the primary protostar. 

In the Z ^ 10 models, fragmentation occurs everywhere in the cloud, and no 

stable rotation-supported disc is formed. Therefore, although the protostars interact with 
each other, they essentially evolve independently. Consequently, several protostars of similar 
mass appear around the centre of the collapsing cloud. On the other hand, the Z > 10 
models yield a single protostar enclosed by a stable rotation-supported disc, which occa¬ 
sionally fragments. The disc develops a spiral structure that effectively transfers the angular 
momentum outwards, suppressing frequent fragmentation and new protostar formation. In 
addition, if a fragment or protostar forms in the disc, it eventually falls onto the central 
protostar. Thus, only the central protostar alone can signihcantly increase its mass in the 
Z > 10“'^ Zq models. 

5 DISCUSSION 

5.1 Star Formation Mode 

The calculation results indicate that the star formation mode considerably differs in the 
lower- and higher-metallicity clouds. In the lower-metallicity clouds (Z ^ 10“^), no stable 
disc is formed, and frequent fragmentation leads to abundant protostars. Although some pro¬ 
tostars are ejected from the central region, the remainder form a small cluster of ~ 10 — 20 
protostars, as shown in Figures |4}l7l The star formation process in the Z = 10~^, 10"^ 


10 ^ Zm models is verv similar to that in the primordial fZ = Oj case ( 

Clark et al. 2008. 

2nila 

b: Smith et al. 2011: 

Greif et al. 2011. 2012: Machida & DoJboiOT). On the other hand. 

the higher-metallicity clouds (Z > 10“'^ Zq models) yield a single protostai 
cumstellar disc, as shown in Figures IHlfTTl Although fragmentation occasioi 

disc, each fragment eventually falls onto the central protostar. The star foi 

the Z = 10“^, 10“^ and 10“^ Zr^ models is similar to the present-dav case ( 

enclosed by a cir- 

aally occurs in the 

'mation process in 

Vorobvov & Basu 

o 

o 

CM 

201ol: Machida et al. 

2011a: Tsnkamoto & Machida 2011: Tsnkamoto et al. 2OI3I1. in 


which fragmentation rarely occurs in the disc and the fragments usually fall onto the central 
protostar. The different star formation modes between clouds of lower-metallicity (small 
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18 M.N. Machida & T. Nakamura 

stellar cluster) and higher-metallicity (star-disc system) are presumably due to differences 
in both the lifetime of the hrst core and mass accretion rate. 


5.1.1 Effects of First Hydrostatic Core 


We hrst discuss the effects of the hrst core on the star formation process. In present-day 
star formation (or the solar-metallicity case), after the collapsing gas becomes optically 
thick to dust thermal emission at n ~ 10^° cm~^ , the gas beh aves adiabatically and the hrst 


(hydrostatic) core forms in the collapsing clond f Larson 


1969 ). The size and mass of the hrst 


core is 0.01 AU and 0.01 Mq, respectivel v flMasunaga fc Inutsukall2000l) . depending on 


the rotation rate of the star-forming cloud (jSaigo fc Tomisakal l2006l) . As the gas accretes 
onto the hrst core, its central density gradually increases and the second collapse, which is 
induced by dissociation of the molecular hydrogen, occurs only in the central part of the hrst 
core at Uc ~ 10^®“^® cm“^. The gas again becomes adiabatic and the protostar (or second 
hydrostatic core) of radius of ~ Rq and mass of ~ 10“^ Mq, hnally appears. The thermal 
evolution of a non-rotating present-day cloud is also reveale d in Figure [H 


Since the rotation of the molecular cloud is non-negligible (iGoodman et ah 


1993 


Caselli et ah 


20021 ). the hrst core is partly supported by the rotation (the thermal pressure is the primary 


support). With the rotation, the hrst core develops a spiral structure that ehectively trans¬ 


fers t 


le angular moment outwards, suppressing fragmentat ion before the protostar forma¬ 


tion flBate 


1998 


Machida et al. 


2008al: 


Tomida et al. 


20131) . After the protostar formation. 


the hrst core remains and evolves into a circumstellar (or rotation-snpported) disc. Only 
the central part of the hrst core alone becomes the protostar, whereas the remainder ro¬ 
tates around the protostar. Thus, the protostar at its formation epoch is alre ady enclosed 


by a l arge-scale disc- 


2011 


Machida et al 


ike structure evolved from the hrst core or its remnant (IBate 


2010 a : 


Machida fc Matsnmoto 


1998, 


2 OI 2 I ). Since the remnant of the hrst 


core or circumstellar disc is more massive than the pro tostar just after the 


tion, fragmentation can occasionally occur in the disc fjinutsuka et al. 


2010 


protostar 


Inutsuka 


orma- 


201211 . 


However, because the fragments usually fall ont o the hrst formed protostar, a single 


star remains in the ear 


Vorobyov fc Basu 


2012 


2010 - 


Tsnkamoto et al. 


y 


main accretion phase f Vorobvov fc Basu 


Machida et al. 


2011 b 


2001 : 


Tsukamoto fc Machida 


proto- 


Walch et al. 


2011 


2009: 


Walch et al 


20131 ). The present-day star formation process has been well investi¬ 


gated, and the above scenario is widely accepted. In addition, some hrst core candidates have 
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been conhrmec 

bv ( 

)bservations ( 

Chen et ah 

2010 

; 

Unoch et ah 

2010; 

Dunham et ah 

2011; 

Belloche et ah 

2011 

Pineda et ah 

2011; 

Chen et ah 

2012 ! 

Pezzuto et ah 

2012; 

Tsitali et ah 

2013; 

Hirano & Liu 

2014 

)• 


On the other hand, the lack of dust prevents the hrst core formation in the primordial 
collapsing cloud. As seen in Figured] the gas temperature gradually increases in the range 
of lO^cm”^ ^ 10^°cm“^ with a polytropic index of 7 ~ 1.1 in the primordial {Z = 0) 
clouds, but suddenly increases at n ~ 10^°cm“^ in the present-day {Z = Ze) clouds. 
Consequently, in the primordial case, no signihcant structure develops during the protostar 
formation epoch (or gas-collapsing phase). In addition, clouds with solar metallicity virtually 
cease collapsing when the hrst core forms. The rotational timescale then becomes shorter 
than the collapsing timescale, and rotational motion dominates in the hrst core. On the 
other hand, in primordial clouds, where the collapsing timescale is always shorter than the 
rotational timescal e with 7 ~ 1 . 1 , the gas continues to fall towards the cloud centre or 
onto the protostar (jSaigo et al.l 120001) . Note that although fragmentation can occur prior to 


the protostar formation if the primordial star-forming cloud has a sufhcien 


momentum, the gas collapse continues in each fragment flMachida et ah 


ly high angular 


20081). Consequently, 


the structure around the protostar during the protostar formation epoch considerably differs 
between the present-day and primordial clouds. The protostar is enclosed by a disc in present- 
day clouds, but is surrounded by radially infalling gas in primordial clouds. The difference in 
circumstellar environment during the protostar formation epoch influences further evolution 
of the protostar and circumstellar region. 



present the cloud evolution prior to the protostar formation in Figure ITSl We observed that 
the hrst core is clearly established in the Z = 10“^ and 10“^ Zq models. At its formation 
epoch, the size of the hrst core is ~ 5 AU in the Z = 10’ model and ~ 3AU in the 
Z = 10“^ Zq model. In both models, the hrst core increases in size as the gas accretes on it, 
and a spiral structure develops. Moreover, fragmentation occurs before the cloud maximum 
density reaches the protostellar density in the Z = 10’ model. Consequently, in both 
Z = 10“^ and 10“^ Zq models, the protostar at its formation epoch is enclosed by a rotating 
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20 M.N. Machida & T. Nakamura 

disc or spiral structure. The sizes of these structures are ~ 20 AU Zq model) and 

~ 5AU (10“^Zq model), respectively. 

As seen in Figure [T51 the gas in the Z = Zq model continues to collapse for tpg ^ 
—5 yr, and the hrst core appears just before the protostar formation at fps — —5 yr. The 
hrst core forms at high gas density of ~ 10^® cm“^, and is as small as ^ 1 AU. The hrst core 
transforms into a ring-like structure that fragments into two protostars. In the Z = 10’ 
model, no clear hrst core appears until the protostar formation. Thus, during the formation 
epoch, the protostar has no disc-like or spiral structure around it. Evidently, protostars 
in clouds of higher metallicity {Z > 10 and low-metallicity {Z ^ 10 -"^g) form 

and evolve in rather diherent circumstellar environments. It appears that the hrst core 
does not signihcantly ahect the main accretion phase in clouds with metallicity of Z ^ 
1 O-"Z0 . According to Figures HlfTTl the star formation process changes between Z = 10 ^ 
and 10““^ Zq. Fragmentation is frequent and leads to many protostars with no stable disc 
formation in the Z ^ 10~‘^ Zq models, whereas the Z > 10“^ Zq models generate a single 
protostar enclosed by a stable rotation-supported disc. This diherence appears to be related 
to the existence of the hrst core. 


5.J.3 Mass Accretion Rate and Disc Stability 

In this subsection, we analytically discuss the relation between the disk stability and mass 
accretion rate. The star formation process is expected to depend also on the mass accretio n 
rates. As reported in previous studies fe.g. iHosokawa fc Omukaill2009bl:iDopcke et al.ll2013[) . 


the mass accretion rate onto the circumstellar region decreases as the cloud metallicity 
increases. This phenomenon is attributed to the cloud temperature. The metal-rich cloud 
has a lower temperature and the mass accretion rate is rou ghly proporti onal to where 
Tci is the temperature of the initial star-forming cloud flLarsonl 1200311 . As shown in ^ 
fragmentation frequently occurs around the cloud centre in the lower-metallicity clouds, 
whereas higher-metallicity clouds rarely fragment. The stability of the circumstellar disc 
evolving around the protostar is considered to be determined by Toomre’s Q parameter 
Q = (csk)/{ nGT.), w here c^, n and S are the sound speed, epicyclic frequency and surface 
density, respectively fjToomrelll9641) . The epicyclic frequency k can be replaced by the Kepler 


angular velocity hlx [= (GMps/r^)^A] -when Kepler rotation is realized in the disc. For 
simplicity, we discuss the disc stability in terms of the parameter 
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ttGS’ 


(4) 


and fragmentation is considered to occur when Q < 1. In our calculation, since some models 
develop no clear rotation-supported disc, we cannot assess the adequacy of Toomre’s analysis. 
We tentatively use it to roughly analyse the stability around the protostar. To simplify the 
problem, we also assume constant sound speed. This assumption is applicable because we 
investigate the evolution of the circumstellar environment over a sufficiently short time 
(~ 100yr after the hrst protostar formation). During this early phase, protostellar heating 
can be ignored beca use the environmental temperature is initially high near the protostar 


fjOmnkai et al 


2010h . 


Equation (j4]) implies that the central star tends to stabilize the disc on a Keplerian 
timescale fxep = 27r/r2K, while the increasing disc surface density tends to destabilize the 
disc on a disc growth timescale tgrow = Mdisc/M. Here, Mdisk and M are the disc mass 
and mass accretion rate onto the central region, respectively. During the main accretion 
phase, the disc growth timescale may be shorter than the Keplerian timescale. In such a 
case, Q continues to decrease and fragmentation is expected to be induced. To very roughly 
investigate the stability of the disc, we compared the Keplerian timescale fxep with the 
disc growth timescale tgrow Figure fTHl plots these timescales in the Z = 10“®, 10“"^, 10“^ 
and 10~^ Zq models for different protostellar masses, namely Mps = 1 (all models), 0.1 
(Z = 10“"^, 10“^ and 10~^ Zq models) and 10 Mq (Z = 10~^ Zq model). The mass accretion 
rate is estimated with a sink (see ^ and Fig. |2]), and the disc mass is assumed as 1% of 
the central protostellar mass, that is, Mdisc = 0.01 Mpg. The protostellar radii Rps at each 
epoch (Fig. [3]) are also plotted in this hgure. Since the region of r < Rps dehnes the inner 
region of the protostar, it is excluded from the disc stability analysis (we call this region the 
forbidden region). 


Outside the protostar (r > Rps), we assume that the disc is stable when txep < t 


grow 


and unstable when txep > igmw In the latter case, it is expected that the disc fragments 
before the central star restores the disc surface density to a stable state. In other words, the 
denominator of equation (j4]) rapidly increases and Q < 1 is soon realized when txep > igrow 
Figure [I6] indicates that when Mps = IMq (red lines), the whole region is gravitationally 
unstable in the Z = 10“® Zq model. In the Z = 10“^ Zq model, a small stable region 
exists in the range of 0.25 AU ^ r ^ 1.5 AU. In the Z = 10“^ and 10’ models, the 
stable region enlarges to 0.05 AU ^ t ^ 10 — 30 AU. Thus, the stable region expands as 
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the metallicity increases. Moreover, in the Z = 10“^ Zq model, the stable region appears 
only when Mps > IMq. This very rongh estimate (there are many ambiguous parameters) 
suggests that fragmentation occurs without stable disc formation in the lower-metallicity 
models because the disc is destabilized by the high mass accretion rate before the disc 
reaches a stable state. On the other hand. Figure [16] indicates that a stable circumstellar 
disc is maintained in the Z > 10“^ Zq models even when the protostellar mass is as small 
as 0.1 Mq (black lines). 

Figure HTj shows the density and velocity distributions at tps — 100 yr in the Z = 10“®, 
10“^, 10“^ and 10“^ Zq models. At this epoch, the most massive protostars are 1.5 Mq 
(Z = 10-5 Zq), IMq (Z = 10-^ Zq), 0.04 Mq (Z = 10-3 Zq) and 0.04 Mq (Z = 10-^ Zq). 
The 10-5 Zq and 10' models develop no stable disc, and many fragments form around 
the cloud centre; in contrast, a disc-like structure is conhrmed in the 10-^ and 10’ 
models. Thus, the simulation results approximately agree with our simple estimate shown in 
Figure [T6l In summary, a high mass accretion rate leads to vigorous fragmentation, whereas 
a lower rate permits a stable rotation-supported disc formation. As the protostellar mass 
increases, a stable disc should form even in the lower-metallicity clouds. However, at least 
during the very early phase of the star formation, fragmentation frequently occurs and many 
protostars form around the centre of the collapsing cloud. Although further calculations are 
necessary to determine the fate of protostars, the protostellar environments of lower- and 
higher-metallicity clouds are considerably different. 
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5.1.3 Gas Temperature and Disk Stability 

Although the Z = 0, 10“® and 10“® Zq models have almost the same accretion rate (Fig. |2]), 
the number of fragments is larger in the Z = 10 model than in the Z = 0 and 10 ® Zq 
models (Fig. [12]). This is expected to be due to the difference in thermal evolution. 

As shown in Figure H] fragmentation tends to occur in a lower-temperature region in 
the Z = 10“® Zq model. On the other hand, in the Z = 10“® Zq model, the cloud centre 
has a relatively low temperature and the entire fragmentation region is embedded in the 


low-temperature region, as shown in Figure [HI The temperature decrea se in the high-c 


region (10^^ cm ^ ^ n ^ 10^® cm ^ of Fig. [T]) is caused by dust cooling (jOmukai et al. 


ensit y 


20101 ). 


Figures [H [4] and [6] indicate that the gas temperature around the protostar (or fragmentation 
region) is lower in the Z = 10“^ Zq model than in the Z = 0 and 10“® Zq models. The lower 
temperature and smaller sound speed decrease the Q parameter and make the disk unstable. 
A temperature gap in the high-density region is also seen in the Z > 10 models. How¬ 
ever, such models have a relatively small mass accretion rate (Fig. [2]). Both a considerably 
large accretion rate and lower-temperature environment are realized in model Z = 10 
Thus, it is expected that the disk in the Z = 10“^ Zq m odel is the most uns t able a nd many 
fragments will appear. A similar result was presented in [Tanaka fc Omukai[ (120141) . 


5.2 Effects of Initial Conditions 


In this study, as the initial state, we adopted the Bonner-Ebert sphere, which is in a near¬ 
equilibrium state. The properties of Bonner-Ebert clouds are determined by the central 
density and temperature, and the cloud temperature is determined by the cloud metallicity 
(see [|2|). Therefore, the cloud masses and sizes differ between models, as listed in Ta ble [H 


Instead, as in previous studies (iDopcke et al. 


2011 


2013 


Safranek-Shrader et al. 


20141 ) , even 


when the star-forming clouds have different metallicities, we can set an identical cloud as 
the initial condition in which the initial cloud has the same mass and size independent of 
its metallicity. However, in such a setup, the cloud stabilities differ between models because 
the initial cloud is in a non-equilibrium state (the equilibrium state is determined by the 
cloud temperature or metallicity). In some previous studies, although the (isothermal) tem¬ 
perature was assumed to be identical in initial clouds with different metallicities, the cloud 
temperature would rapidly change and each cloud would have different thermal energies (or 
different stability) just after the cloud began to collapse. 
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The (initial) cloud stability, defined as the ratio of thermal to gravitational energy ao, 
significantly affects the fragmentation condition in the colla psing cloud. Clouds with low er 


thermal energies (smaller pq) are prone to fragmentation flTsuribe Sz Inutsuka 


1999a 


M). 


Thus, it is expected that under identical initial conditions, a higher-metallicity cloud will 
more likely fragment because the thermal energy is reduced by efficient cooling. In this 
scenario, the initial condition determines the outcome (or frequency of fragmentation). Es¬ 
pecially, unless the high-density gas region or the first core is not resolved with a sink, 
frequent fragmentation is expected in higher-metallicity clouds because Pq is initially small, 
as seen in previous studies. 

It is very difficult to determine adequate initial conditions for the general star formation 
process. The initial conditions of a prestellar cloud wi th zero metallicity may be determined 
by cosmological simulations (e.g. 


Hirano e 


ah 


201411. where as those of a solar-metallicity 


Motteetah 


19981 ). In both cases, i nitial clouds 


cloud are deduced from observations (e.g. 
are c onsidered to be well descri bed by the Bonne r-Ebert density profile (see 


Abel et ah 


2002 


for the primordial case and 


Alves et al. 


12001 


for the solar metallicity case). However, 


the properties of prestellar clouds (or adequate initial conditions) with 0 < Z < Zq cannot 
be determined. Although it seems natural to assume that such clouds also exist in a near- 
equilibrium state with a Bonnor-Ebert density profile, their properties cannot be confirmed. 

The fragmentation condition depends not only on the thermal energy but also on the ini¬ 
tial cloud rotation. Clouds with initially higher rotational energy (or larger f3o) are more likely 


to fragment. In our previous studies, we in vestigated the effects q 


star formation pr ocess in both primordial flMachida et al. 


solar -metallicity flMachida et al. 


2005b; 


Machida et al. 


2008c: 


2008 


the initial rota 


M 


2010 a : 


achida fc Doi 


Ion on the 


20131 and 


Machida Sz Matsumoto 


2 OI 2 I ) clouds. These studies showed that, independent of (do, solar-metallicity clouds develop 


a stable rotation-supported disc, whereas primordial clouds undergo vigorous fragmentation, 
leading to numerous protostars with no stable discs. Identical results were obtained in the 
primordial and solar-metallicity clouds investigated in the present study. To better under¬ 
stand the accretion phase of star formation in clouds with different metallicities, we may 
need to investigate the cloud evolution with different (3o. However, such a parameter survey 
is currently precluded by the very long calculation time required to investigate the accre¬ 
tion phase resolving the protostar. The typical computational time of a single model in this 


study was six months. Although we expect that changing (3o w ill not qualitatively c 


range 


the cloud evolution, as previously found in primordial clouds flMachida Sz Doil l2013l) and 
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solar-metallicity clouds with a sink flMachida et al.ll2010al) . the effect of rotation in clouds 


with different metallicities requires further investigation. 

In this study, we limited the number of models and focused on the effect of thermal 
evolution, which is dependent on cloud metallicity. Therefore, we constrained the initial 
cloud to a near-equilibrium condition and a constant ratio of rotational to gravitational 
energy. However, determining the correct initial conditions is essential for understanding the 
evolution of clouds with different metallicities. Thus, in future work, we must investigate 
cloud evolution over a large parameter space of ao and /Sq. 


5.3 Effects of Magnetic Field 


In this study, we ignored the effects of the magnetic field. However, the magnetic field 
has a large impact on star formation, because it effectively transfer s the angular momen¬ 


tum, thereby suppressing cl oud fragmentation if sufficiently strong flMachida et ah 


Hennebelle fc Tevssier 


2008a; 


20081) . In addition, a po rtion of the infa 


tostellar outflow driven by the Lorentz force flTornisaka 


2002 


l ing gas is ejected by t 


formation is also suppressed by magnetic braking flMellon fc Li 


Machida et al 


200 


2009; 


le pro- 


2008b). Disc 


Machida fc Doi 


2013|). 


To investigate the effects of the magnetic field in clouds with different metallicities, 
we require both the magnetic field strength of the star-forming cloud and the dissipation 
process of the magnetic field in each cloud (or each metallicity). The magnetic field strength 
of solar-metallic ity clouds is determined from observations of nearby star-forming regions 


fjCrntcherl 1 19991 ). On the other hand, the field strengths of clouds with Z Zq are less 


easily determined because there is no observational evidence. 

In addition, the magnetic field can dissipate and weaken in a collapsing cloud. In clouds 


Nakano et al. 

2002 ; 

Machida et al. 

2007) 


do not lose their magnetic fields by dissipation flMaki fc Susa 


le other hand, col 


2004 


apsing primordial clouds 


20071) . The difference in 


the magnetic dissipation process between solar-metallicity and primordial clouds is primar¬ 
ily due to the different ionization degree. Clouds with solar metallicity are relatively poorly 
ionized because of their lower temperature and abundant dust. In contrast, a primordial gas 
has a high temperature and lack of dust, and is ionized to a relatively high degree. Since 
the temperature of the collapsing cloud and dust abundance depend on the metallicity, the 
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dissipation degree of the magnetic field should also be affected by the metallicity. Although 
the initial strength and dissipation degree of the magnetic held are not easily determined, in¬ 
vestigating the effects of the magnetic held on star formation processes in clouds of diherent 
metallicities is an important future task. 


5.4 Effect of Stellar Feedback 


In this study, we ignored the ehect of stellar feedback that can suppress disk fragmentation 
through heating of the circumstellar disk by the protostellar luminosity and stabilization 
against gravitational instability. In this subsection, we discuss the ehect of protostellar feed¬ 
back on fragmentation. 

In present-day low-mass star formation, the prot ostellar feedback ahects the disk evo¬ 


lution and fragmentation process. lOhner et ah 


( 20091 ) 


investigated low-mass star formation 


in a turbulent cl oud and sh o wed that protostellar 


tation (see also 


Bate 


2009 


Krnmholz et al 


2010 ) 


e edback tends to su p press disk fragmen 


Stamatellos et al. 


( 20121) showed that 


episodic accretion alleviates the suppression of disk fragmentation due to pro tostellar feed 


back because fragmentation can occur in a quiescent phase between outbursts (iLomax et al 


2014J ). Thus, in the present-day star formation process, although protostellar feedback cannot 
completely suppress disk fragmentation, it plays an important role for the disk evolution. 


On the other hand, t he ehect of s 
lower-metallicity clouds. IStacv et al. 


ellar fe edback is not very signihcant in primordial and 


(120121) showed that radiative feed back does n o 


vent disk fragmentation and multiple stars appear in primordial clouds. ISmith et al, 


pre 


(1201 ih 


calculated the evolution of primordial minihalos including the luminosity from accreting pro¬ 
tostars, and showed that accretion luminosity does not prevent disk fragmentation. This is 
because the heating from the accretion luminosity contributes little to the disk thermal evo¬ 


lution; the compression heating dominates other heating pr ocesses. To estimate 


protostellar feedback in clouds with diherent metallicities. 


Omukai et al 


he ehect of 


(120101) compared 


the cloud gas temperature before heating (or before the protostar formation) with that 
heated by the protostar, and showed that the former dominates the latter in clouds with 
Z ^ 10“^ Zq. In lower-metallicity clouds, since the gas temperature is intrinsically high due 
to the dehcit of efhcient coolant, protostellar feedback does not play a signihcant role for 
disk evolution and fragmentation. On the other hand, in higher-metallicity clouds, since the 
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gas temperature is relatively low, the protostellar luminosity effectively heats the disk and 
suppresses fragmentation. 

In our study, fragmentation frequently occurs in lower-metallicity clouds, while it rarely 
occurs in higher-metallicity clouds. Since previous studies indicate that the effect of proto¬ 
stellar feedback is not significant in lower-metallicity clouds, we expect t hat frequent frag¬ 


mentation occurs in lo wer-metallicity clouds, as seen in primordial clouds flGreif et al.l 


2012 


Stacy fc Rromm 


20141) . In addition, at the end of the calculation, the protostars are found 


to have masses of ~ 1 Mq in lower-metallicity clouds {Z ^ 10’ -"^ 0 ) and < 0.1 Mq in higher- 
metallicity clouds {Z > 10~'^Zq). Thus, it is expected that the protostars exhibit minimal 
radiative effects during the calculation. However, to quantitatively understand fragmenta¬ 
tion properties such as stellar mass and the number of fragments, we need to calculate the 
cloud evolution with stellar feedback. 


5.5 Can Ejected Low-metal Stars Be Observed? 

In clouds with Z ^ 10“^ Zq, many protostars form by fragmentation and low-mass protostars 
are preferentially ejected from the cloud centre. In addition, low-mass protostars with masses 
of <C 1 Mq exist even around the centre of the cloud at the end of the calculation. These 
protostars may evolve into zero- (or extremely low-) metallicity {Z ^ 10 -"^ 0 ) and low- 
mass (M ^ 0.8 Mq) stars that can survive until the present day and may be observable. In 
this subsection, we discuss the observation probability of such zero-metallicity or extremely 
metal poor stars in our galaxy. 

We evolved each cloud ~ 100 yr after the first protostar formation and found that frag¬ 
mentation frequently occurs and many protostars appear in zero- and lower-metallicity en¬ 
vironments. Some protostars have a mass of <C 1 Mq at the end of the calculation. However, 
it is expected that such low-mass protostars increase their mass to reach > 1 Mq in a further 
evolution stage because the zero- (or low-) metallicity host cloud has a sufficient mass to 
sustain the mass accretion for a long duration. In addition, the mass accretion rate is con¬ 
siderably high in zero- and lower-metallicity clouds. Thus, we expect that it is considerably 
difficult for a protostar around the cloud centre to keep its mass within ^ 0.8 Mq. The stars 
with a mass of > 0.8 Mq would have died. 

On the other hand, since the ejected protostar moves away from the cloud centre and mass 
accretion is expected to gradually weaken, such stars may have a mass of ^ 0.8 Mq. Although 
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we need to calculate the cloud in further evolution stages to determine the hnal stellar 
mass and the number of zero-metallicity (or extremely metal poor) stars still surviving, 
we probatively estimate the number of observable zero-metallicity stars in our galaxy. We 
dehne the number of low-mass stars with a mass of ^ 0.8 Mq in a single minihalo as A^mh- 
The minihalo is assumed to have a baryonic mass of Mmh = 10® Mq, while the mass of 
our galaxy is Mgai ~ 10^^ M©. The galaxy evolves by the merger of minihalos and we 
simply assume that about 10® minihalos are gathered to form our galaxy. Thus, our galaxy 
contains N x 10® zero-metallicity stars. Therefore, the number ratio of zero-metallicity to 
Population I and II stars is iV x 10“®. Simply assuming iV = 1 implies that we can hnd 
a zero-metallicity star randomly in the observation of one million stars in our galaxy. In 
other words, we need to observe about 10® stars in our galaxy to hnd a zero-metallicity 
star. Although sev eral hundred thousand stars have been spectroscopically identihed (e.g. 


Yannv et al. 


2009f) . further observations may be necessary to hnd zero-metallicity stars. 


Although the observational probability of zero-metallicity stars strongly depends on N, 
we expect that N is not large. Since the ejected protostar can acquire its mass from the 
infalling envelope even after being ejected from the central region, it is difficult to keep 
its mass in the range of < 0.8 Mq. In addition, some low-mass protostars are expected 
to be further ejected from minihalos and may drift in intergalactic sp ace. Moreoyer, an 


initially zero-metallicity sta r can be pollu 
or interstellar medium (e.g. 


Komiva et al. 


ed by a binary comp anion (e.g. 


Suda et al. 


2010 


Johnson 


2004h 


20141) . and may be observed as an 


extremely metal poor star. Although a more detailed study is necessary to determine the 
observational probability of zero- (and extremely low-) metallicity stars, we may observe 
them in the near future. 


6 SUMMARY 

In this study, we investigated the early main accretion phase of star formation in clouds with 
diherent metallicities in the range of 0 ^ ^ Zq. We calculated the cloud evolution by a 

three-dimensional nested grid code, in which the thermal evolution given by the one-zone 
model was used. As the initial state, we adopted the Bonner-Ebert cloud, which exists in 
a near-equilibrium state. Each initial cloud had different metallicities and temperatures. In 
other words, the initial clouds differed in mass and size, but were assigned the same ratio of 
thermal (cto) and rotational (/Sq) energy to gravitational energy. Starting from the prestellar 
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stage, we calculated the cloud evolution for ~ 100 yr after the first protostar formation, 
resolving the protostellar radius without introducing a sink. To adequately calculate the 
circumstellar environment, we adopted a protostellar model, which relates the protostellar 
radius to the protostellar mass. 

We showed that the star formation process considerably differs between the lower- and 
higher-metallicity clouds. In the lower-metallicity clouds {Z ^ 10“'^ Zq), fragmentation fre¬ 
quently occurs without a stable disc formation and over ~ 10 protostars appear at the 
same epoch around the cloud centre. The number of protostars is slightly smaller in the 
Z = 10“'^ Zq model than in the Z < 10“'^ Zq models. Some of the protostars fall onto a 
more massive protostar. Thus, a massive protostar acquires its mass by two processes: gas 
accretion and the merger of less massive protostars. The most massive protostar has a mass 
of ~ 1 Mq at the end of the calculation. On the other hand, less massive protostars with 
masses of ~ 0.01 Mq are preferentially ejected from the central region by mutual gravita¬ 
tional interactions. They are expected to become lower-metallicity low-mass stars. At the 
end of the calculation, the cloud centre is surrounded by 7-20 protostars forming a small 
stellar cluster. 

On the other hand, in the higher-metallicity clouds of Z > 10-4 Zq, fragmentation 
rarely occurs and a single protostar appears. Immediately after its formation, the protostar is 
enclosed by a circumstellar disc, which is supported by rotation. Since the disc and protostar 
have similar masses, the disc fragments and develops a few clumps. The clumps fall onto 
the central protostar before evolving into separate protostars. The disc develops spiral arms, 
which effectively transfer the angular momentum outwards. At the end of the calculation, in 
these models, a single protostar is enclosed by a circumstellar disc, which evokes the classical 
star formation process. 

In summary, the Z ^ 10 models yields a small stellar cluster with no stable disc, 
whereas Z > IQ-^Zq produces a protostar surrounded by a circumstellar disc. In other 
words, lower-metallicity clouds fragment everywhere, and protostars independently evolve 
with mutual interaction but no stable disc formation. Their higher-metallicity counterparts 
yield a single protostar enclosed by a large-scale circumstellar disc, which only occasionally 
fragments. This evolutionary difference is caused by metallicity-dependent thermal evolu¬ 
tion and the mass accretion rate. In the Z > Zq models, the protostar formation is 
proceeded by a clear hrst core formation, and the hrst core ultimately becomes the circum¬ 
stellar disc. The first core (or circumstellar disc) stabilizes the circumstellar environment 
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against fragmentation with effective outward transfer of angular momentum. This scenario 
is consistent with the present-day (or solar metallicity) star formation process. On the other 
hand, no (clear) hrst core forms in the Z ^ 10’ models. Thus, the gas directly falls onto 
or near the protostar, and fragmentation frequently occurs without an effective mechanism 
of angular momentum transfer. Therefore, a stellar cluster composed of ~ 10 protostars ap¬ 
pears around the cloud centre. This scenario mirrors the primordial star formation process. 
Furthermore, the high mass accretion rate prevents stable disc formation in the Z ^ 10““^ Zq 
models. The mass accretion rate is moderated in the Z > 10 models, permitting the 
formation of a stable circumstellar disc. 

In this study, we showed that the star formation process transition occurs at Z = 10“'^ — 
10“^ Zq. However, we calculated the cloud evolution over a very limited time: ~ 100 yr after 
the hrst protostar formation. We also ignored the magnetic held, which is expected to have 
a large impact on star formation. To properly clarify the ehect of metal abundance on the 
star formation process, we must calculate the cloud evolution over a longer time and include 
the magnetic held and its dissipation process by imposing adequate initial conditions. 
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Table 1. Models simulated in this study and their results. Column 1 gives the model name. Column 2 gives the cloud metallicity 
{Z). Columns 3-5 give the temperature (Td), mass (Md) and radius (Rd) of the initial cloud. Column 6 gives the protostellar 
density (rips). Column 7 gives the cell width of the finest grid. 


Model 

Z 

Tci [K] 

Mel [Me] 

Rcl [AU] 

rips [cm 3] 

h [AU] 

1 

0 

197 

2500 

3.8 X lO® 

3.4 X lOil" 

3.8 X 10-2 

2 

10-6 Zq 

195 

2460 

3.8 X lO® 

3.8 X lOil 

3.8 X 10-2 

3 

10-5 Zq 

190 

2370 

3.7 X lO® 

4.3 X lOil 

3.9 X 10-2 

4 

10-4 Zq 

154 

1730 

3.4 X 10® 

3.4 X lOi® 

8.6 X 10-3 

5 

10-3 Zq 

34 

180 

1.6 X lO® 

8.5 X lOi® 

4.3 X 10-3 

6 

10-2 Zq 

18 

53 

1.0 X 10® 

9.1 X lOi® 

2.2 X 10-3 

7 

10-1 Zq 

20 

62 

1.1 X 10® 

1.7 X 1020 

1.1 X 10-3 

8 

Zq 

11 

9.5 

4.8 X 104 

2.2 X 1020 

1.1 X 10-3 
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Figure 1. Temperature of collapsing gas clouds with different metallicities against the number density. The metallicity of each 
thermal evolution is indicated in the figure. Two evolutionary tracks are plotted for each metallicity: one given by the one-zone 
calculation (black) and the other by the one-zone result added to the protostellar model in a high-density region (red). 



Figure 2. Mass accretion rate onto a sink versus protostellar mass. Results are plotted for different metallicities (indicated 
at the right of the plots). 
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1/3 

Figure 3. Protostellar radius versus protostellar mass for different meta llicities. The r e lation Rps oc Mpg is also plotted. 
Protostellar radius for Z = 10“^, 10“"^, 10“^, 10~^ and 1 Zq calculated in lOmukai et all i|2010l ') are also plotted as thin lines. 
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Figure 4. Density (upper panels) and temperature (lower panels) distributions on the equatorial plane in two models: Z = 0 
(left) and 10“® Zq (right). The elapsed time after the calculation starts t and that after the protostar formation tps are 
described in each upper panel. 
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Figure 5. Time sequence of density distribution on the equatorial plane in the Z = 10“^ Zq model. The elapsed time after 
the cloud begins to collapse t is described in panel (a) and that after the protostar formation tps is described in each panel. 
The box size is different in each panel. 
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Figure 6. Temperature distribution in Z = 10 ® (left) and 10 Zq (right). The elapsed time after the protostar formation 
ips is described in each panel. The epochs of the left and right panels correspond to Fig. (5)1 and Fig.[7lf, respectively. 
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Figure 7. Same as Fig. [5] for the Z = 10 ^ Zq model. 
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Figure 10. Same as Fig. [5]for the Z = 10 ^ Zq model. 
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Figure 11. Density distribution on the equatorial plane in the Z = 1 Zq model. The elapsed time after the calculation starts 
t and that after the protostar formation tps are described. 
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Figure 12. Number of protostars versus the elapsed time after the first protostar formation. Results are plotted for all models. 
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Figure 13. Mass of the most massive protostar versus the elapsed time after the first protostar formation. Results are plotted 
for all models. 
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Figure 14. Mass of the first, second, fourth, sixth, eighth and tenth most massive protostars in models Z = 0, 10 10 

10“^, 10-3 and lO-^ Zq versus the elapsed time after the first protostar formation. The total mass of protostars is also plotted 
by a broken line in each panel. 
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Figure 15. Time sequence (rows) of cloud evolution on the equatorial plane in the Z = 10~^, 10“^, 10”'^ and 10~^ Zq 
models. The density distribution (colour) and velocity vectors (arrows) on the equatorial plane are plotted in each panel. The 
elapsed time after the cloud begins to collapse t and that after the protostar formation tps are also described in each panel. 
The box size is different in each panel. 
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Figure 16. Disc growth tgrow and Keplerian timescale t^ep for -^ps = 1 (red), 0.1 (black) and 10 Mq (blue) versus the radius 
from the protostar in the Z = 10~^, 10“^, 10“^ and 10“^ Zq models. The protostellar radius Rps in each epoch is also plotted 
in each panel. The unstable and stable regions are defined as tgrow < ^Kep and tgrow > ^Kep outside the protostellar radius 
r > -Rps, respectively. 
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Figure 17. Density (colour) and velocity (arrow) distributions on the equatorial plane at tps — 100yr in the Z = 10“^, 10“^, 
10“^ and 10“^ Zq models. The elapsed time after the cloud begins to collapse t and that after the first protostar formation 
£ps are described in each panel. The spatial scale is different in each panel. 
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APPENDIX A: TEST CALCULATION FOR SPHERICAL CLOUDS WITH 
PROTOSTELLAR MODEL 


To verify our protostellar model, we calculated a non-rotating cloud with different metal- 
licities as described in ^ Figure lATl shows the density (upper panels) and velocity (lower 
panels) prohles at different epochs in clouds with Z = 10“^ Zq (left panels) and 10“^ Zq 
(right panels). The initial state for the models is a non-rotating cloud. The protostellar sur¬ 
face is marked by a sudden increase in density in a high-density region. The hgure indicates 
that the protostar gradually enlarges over time. In addition, the prohles correspond well 
with the analytical solutions. In the gas-collapsing phase prior to the protostar formation, 
the density is proportional to p oc in the outer region, an d the radial velocity (vr) grad¬ 


ually decreases with decreasing radius near the cloud centre flLarson 


1969 


Omukai fc Nishi 


19981 ). On the other hand, during the gas accretion phase after protostar formation, the 


velocity continuously decreases with decreasing radius and suddenly becomes n,, — 0 at the 
protostellar surface. The density prohle has p oc near the protostar and p oc far 
from the protostar. IWhitworth fc Summers! (119851 ) showed that the dens ity profile around 


the protostar chang es from p oar ^ to p oc r after protostar formation 


1977 


Hunter 


f Larso J 


1969 


Shu 


19771) . These features are clearly seen in Figure |AT] (left panels). 


Immediately before and after the protostar formation, the outer envelope is more com¬ 
plicated in the Z = 10“^ Zq model than in the Z = 10“^ Zq model. The red lines in the 
right panels of Figure lATl show that, in the Z = 10 model, there are two shock fronts 
at r 5Rq ( 0.02 AU) and ~ 60 Rq (~ 0.3 AU). The inner shock corresponds to the pro¬ 

tostar, whereas the outer shock corresponds to the hrst core or its remnant. In a collapsing 
cloud containing abundant metals and dust, after the dust cooling becomes optically thick 
and adiabatic heating dominates the dust cooling, the gas behaves adiabatically in the range 
of 10^° cm“^ ^ 10^® cm“^ (Fig. [p. During the adiabatic contraction phase , the hrst adi¬ 


abatic core (or the hrst core) forms flLarsonlll969c iMasunaga fc Inutsuka 


20001) . However, as 


seen in the right panels of Figure lAll the hrst core disappears several years after the pro¬ 
tostar formation, when the star-forming clond has no angular momentum. Note that when 
the star-forming cloud has angular momentum and contains a greater abundance of metals 
le hrs 


disc 

Bate 

1996; 

Machida et ah 

2010a: 

Bate 

2011: 

Bate et al. 

2014; 

Machida & Matsnmoto 

2011; 

Tomida et ah 

2013; 

Tsukamoto et al. 

201 3I1. 

Thus, the hrst core plays a critical role in 
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the star formation p rocess. In the Z = 10' model, the hrst core exists for a very short 
period (see Fig. 3 in iMachida et al.ll2009bl) and immediately disappears after the protostar 


formation. The effect of the hrst core and cloud metallicity on the star formation is discussed 
in 1 15.1.11 

Figure |A2] shows the time sequence of density and velocity distributions on the equatorial 
plane after the protostar formation in the Z = 10“^ Zq (upper panels) and 10“^ Zq (lower 
panels) models. The central white region bounded by the shock front indicates the protostar. 
The hgure indicates that the protostar gradually enlarges with increasing protostellar mass, 
as seen in Figure IXTl In addition, the surrounding gas density gradually decreases, whereas 
the (negative) radial velocity increases with time. Figures lATl and lA^ land Figs. 4 and 5 in 


Machida fc Doi 


20131) indicate that in the absence of angular momentum, the protostellar 


environment in the proximity of the protostar shows a similar evolution among models with 
different metallicities although the protostellar sizes (and mass accretion rates) considerably 
differ. On the other hand, even in non-rotating clouds there is a great difference in the large 
scale structure 0.1 — 1AU), because the hrst core which is formed in a cloud with higher 
metallicity with a size of ~ 0.1 — lAU surrounds the protostar and greatly changes the 
circumstellar environment in the region of 0.1 — lAU (for details, see 1 15.1.ip . When the 
cloud has angular momentum, the star formation process is largely governed by the presence 
of the hrst core and the mass accretion rate. 

In this appendix, we have shown that our model approximately reproduces protostel¬ 
lar evolution. In ^ we examine the evolution of initially rotating clouds having diherent 
metallicities with the protostellar model. Note that our modelling ignores the protostellar 


feedback ehects of the ambient and infalling gas. However, 


Omukai et al. 


( 20101) pointed out 


that when the protostar is not very massive (^ 10 — 100 Mq), stellar feedback to the accre¬ 
tion how is not signihcant yet. This is because the high-density gas region, which is the focus 
of this study, has a rel atively high temperature in the early star formation stage (see Fig. 21 


m 


Omukai et al. 


20101) . Also, the ehect of the radiation pressure can be ignored in the very 


early phase of the stm forrn ation because the protostellar luminosity is sufficiently small 
(iHosokawa fc Omukaill2009bl) . Note also that since the mass accretion rate is expected to be 
diherent in each protostar after fragmentation, the mass-radius relation shown in Figure [3] 
may not be correctly reproduced in each fragment. In this study, we calculated the cloud 
evolution for only ~ 100 yr after the protostar formation, where the protostars have a mass 
of < 10 Mq. Although more detailed modelling may be necessary to correctly investigate the 
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Figure Al. Density (upper panels) and velocity (lower panels) profiles at five different epochs for the models Z = 10 ^ Zq 
(left panels) and 10“^ Zq (right panels) versus the radius from the cloud centre. 

protostar formation and its effects, we believe that our model can qualitatively investigate 
the circumstellar disc formation and fragmentation process during early star formation. 
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Figure A2. Time sequence of density and velocity distributions on the equatorial plane in two models: Z = 10“^ (upper 
panels) and 10“^ Zq (lower panels). The spatial scales differ between the upper and lower panels. The protostellar mass is 
described in each panel. 
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